#get the command line arguements
args=(commandArgs(TRUE))

#evaluate the arguments
for(i in 1:length(args)) {
	 eval(parse(text=args[[i]]))
}

#set the working directory
setwd("/home1/30/jc227089/Toad_point_process_model/DONE_NAT")

#####-------------------Functions-----------------#####
#The kernel
dcncross<-function(x, u) {  #stuart's cauchy-normal distribution in 2D with v->Inf
	exp(-x^2/(2*u^2))/(2*pi*u^2)
} 

# outputs data.  With or without a plot as well.
output<-function(pop, gen, id, plot=T){
	dname<-paste(id, "pop", gen, ".txt", sep="")
	write.table(pop, dname, sep="\t", row.names=F)
	if (plot==T){
		pname<-paste(id, "plot", gen, ".png", sep="")
		plotter(pop, pname)
	}	
}

# spreads the population over gen generations
spread<-function(pop, gens, pairs, K){ #pairs is a list from pdist.fast
	for (i in 1:gens){
		#if (i%%5==0) output(pop, i, K)
		occp<-subset(pop, pop[,"Pres"]==1 & pop[,"age"]<6) #collect occupied sites
		if (length(occp[,1])<1) {
			print(paste("No more vacant opportunities at generation", i))
			break()	
		}
		potl<-pairs[occp[,"ID"]] #collect relevant parts of pair list
		potl<-do.call("rbind", potl)
		src.ID<-rep(occp[,"ID"], times=occp[,"n.pairs"])
		potl<-cbind(src.ID, potl)
		U<-rep(occp[,"u"], times=occp[,"n.pairs"]) #expand source specific kernel parameters
		#V<-rep(occp[,"v"], times=occp[,"n.pairs"])
		recruits<-dcncross(potl[,"dists"], U) #calculate probabilities of arrival for each pair
		n.inds<-sum(recruits)*K #Calculates number of individuals to draw
		recruits<-tapply(recruits, potl[,"snk.ID"], sum) #vector of probabilities for the multinomial
		recruit.ID<-as.integer(row.names(recruits))
		recruits<-rmultinom(1, n.inds, recruits)
		recruits<-cbind(recruit.ID, recruits)
		colnames(recruits)[2]<-"recruits"
		recruits<-subset(recruits, recruits[,"recruits"]>2)
		pop[match(recruits[,"recruit.ID"], pop[,"ID"]), "Pres"]<-1
		pop[which(pop[,"Pres"]==1), "age"]<-1+pop[which(pop[,"Pres"]==1), "age"]
	}
	output(pop, i, K)	
}

#plots opportunities and colonised populations
plotter<-function(popmatrix, file.name){
	png(filename=file.name, width=7, height=7, units="cm", res=600, pointsize=6)
	plot(popmatrix[,2], popmatrix[,3], xlab="False easting (kms)", ylab="False northing (kms)", pch=19)
	occp<-subset(popmatrix, popmatrix[,"Pres"]==1)
	points(occp[,2], occp[,3], pch=19, col="red")
	dev.off()
}


##Run model##
load("Pairs.RData")
load("pop.tabl.seeded.RData")

spread(tabl,72,pairs,PseudoK)
